/*
**  File: AquaRange.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann
**  VERSION   : 1.2.1
**  DATE      : Mar 27, 1997
**
**  This file contains the Aqua interface to the Range routines.
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"

/****************  function prototypes  ****************/

#include "AquaFuncts.h"
#include "AquaStrucset.h"
#include "MenuFuncts.h"
#include "Functs.h"

/* BEGIN_EXTH */
#include "Range.h"
/*END_EXTH */

#include "AquaRange.h"

/****************  globals  ****************/

Range           *range;                 /* export */



/* BEGIN_EXTH */

/****************  ReadRanges  ****************/
/*
**  Read a Range, and return it.
*/

Range *ReadRanges( Filnam filnam )
{
    char    *fct_name = "ReadRanges";
    Range   *range;
    int     ifail, both_ends;
    
    range = CREATE_NEW( Range, 1);
    if ( !range )
	NOMEM( "range" );
    get_ranges( filnam, range->first_model, range->last_model,
                range->first_residue, range->last_residue,
                &(range->nmodel_ranges), &(range->nresidue_ranges), 
		&both_ends, &ifail );
    if ( both_ends )
	ModOption( opt_rest_sel, "both" );
    else
	ModOption( opt_rest_sel, "one" );

    return( range );
}

/****************  InResidueRange  ****************/
/*
**  Check whether a residue is in 'range' (global variable).
*/

Boolean InResidueRange( Residue *res )
{
    char   res_num[6];

    sprintf( res_num, "%5i", res->res_num_orig );
    return ( (Boolean)check_residue_ranges( res_num, res->chain[0], 
                                      range->first_residue, range->last_residue, 
		                      range->nresidue_ranges ) );
}

/****************  InModelRange  ****************/
/*
**  Check whether a model is in 'range' (global variable).
*/

Boolean InModelRange( const int current_model )
{
    return ( (Boolean)check_model_ranges( current_model, range->first_model, 
                                       range->last_model, range->nmodel_ranges ) );
}

/****************  RangeStrucset  ****************/
/*
**  Apply a Range to a Strucset, i.e. select residues and/or models. 
**  The operation is similar to StrucsetCopy.
**  NOTE: in RangeDistLists and ProcessDistvalFile modelnumbers are stored.....
**  these should really be part of a Strucset...........................
**  and the model-part should be moved from RangeDistLists to here..........
**  The new Strucset should be Init-ialized, but not Create-d.
**  The Range is a global variable, used in routines InResidueRange. 
**  <atom_map> returns a pointer to an array, containing the new sequential
**  atom indices as a function of the old ones.
**  <residue_map> returns a pointer to an array, containing the new sequential
**  residue indices (not the numbers!) as a function of the old ones.
*/

void RangeStrucset( Strucset *new_variable, Strucset *old, int **atom_map, int **residue_map )
{
    char     *fct_name = "RangeStrucset";
    Mol      *old_molptr, *new_molptr;
    Struc    *old_strucptr, *new_strucptr, *nextptr;
    Atom     *old_atomptr;
/*    Atom     *new_atomptr; */
    Coord    *coordptr;
    Residue  *old_resiptr, *new_resiptr;
    int      *atommap, *resimap, *mapptr;
    int      idx, old_count, new_count;

    STATUS( "entering" );

    if ( !StrucsetOK( old, complete_set ) )
        exit( 1 );
    if ( !EMPTYSET( *new_variable ) )
    {
        printf( " * Will not copy into non-empty strucset <%s>\n", new_variable->name );
        exit( 1 );
    }
/*
    if ( old->struc_count )
        FATAL( "RangeStrucset does not work on coordinates (yet)" );
*/
    if ( !old->onemap )
	ERROR( "IGNORING MODEL DEPENDENT ATOM NAME MAPS" );

    old_molptr = old->molecule;

/* allocate the map arrays */

    atommap = CREATE_NEW( int, old->atom_count );
    if ( !atommap )
	NOMEM( "atommap" );
    resimap = CREATE_NEW( int, old_molptr->res_count );
    if ( !resimap )
	NOMEM( "resimap" );

/* count the number of atoms to be selected */

    old_atomptr = old_molptr->atoms;
    old_resiptr = old_molptr->residues;
    old_count = old->atom_count;
    new_count = 0;
    while ( old_count-- )
    {
	idx = FindResidueIndexFromNum( old_atomptr->res_num_new, old_molptr );
	if ( idx < 0 )
	{
	    ERROR( "Index out of range - Code error" );
	    exit( 1 );
	}
	if ( InResidueRange( old_resiptr+idx ) )
	    new_count++;
	old_atomptr++;
    }

/* allocate the necessary memory for the new strucset */

    StrucsetCreate( new_variable, new_count, old->struc_count );
    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );
    
    new_variable->type = sel_str_data;
    new_molptr = new_variable->molecule;

/* copy the Mol (and save the atoms pointer ) */

    MolCopy( new_molptr, old_molptr );

/* copy the Struc's (and save the coords and next_struc pointers ) */

    new_strucptr = new_variable->fst_struc;
    old_strucptr = old->fst_struc;
    while ( new_strucptr ) 
    {
        coordptr = new_strucptr->coords;                        /* save ptrs */
        nextptr = new_strucptr->next_struc;
        memcpy( new_strucptr, old_strucptr, sizeof( Struc ) );  /* copy */
        new_strucptr->coords = coordptr;                        /* restore ptrs */
        new_strucptr->next_struc = nextptr;
        new_strucptr->atom_count = 0;                           /* reset atom count */
        new_strucptr = nextptr;                                 /* goto next structure */
        old_strucptr = old_strucptr->next_struc;
    }

    if ( !StrucsetOK( new_variable, created_set ) )
        exit( 1 );

/* copy the Atom's and set the atom map */

    mapptr = atommap;
    old_atomptr = old_molptr->atoms;
    old_resiptr = old_molptr->residues;
    old_count = old->atom_count;
    new_count = 0;
    while ( old_count-- )
    {
	idx = FindResidueIndexFromNum( old_atomptr->res_num_new, old_molptr );
	if ( idx < 0 )
	{
	    ERROR( "Index out of range - Code error" );
	    exit( 1 );
	}
	if ( InResidueRange( old_resiptr+idx ) )
	{
	    AddAtomToSet( new_variable, old_atomptr );
	    *mapptr++ = new_count++;            /* new sequential atom index */
	}
	else
	    *mapptr++ = -1;
	old_atomptr++;
    }

/* copy the Coord's, using the map just made */

    new_strucptr = new_variable->fst_struc;
    old_strucptr = old->fst_struc;
    while ( new_strucptr ) 
    {
	mapptr = atommap;
	old_count = old->atom_count;
	coordptr = old_strucptr->coords;
	while ( old_count-- )
	{
	    if ( *mapptr++ != -1 )
		AddCoordToSet( new_variable, new_strucptr, coordptr );
	    coordptr++;
	}
        new_strucptr = new_strucptr->next_struc;
        old_strucptr = old_strucptr->next_struc;
    }

/* check the strucset */

    if ( !StrucsetOK( new_variable, complete_set ) )
        exit( 1 );

/* make the residue list and set the residue map */

    MakeResidues( new_variable, true );
    StrucsetStats( new_variable );

    mapptr = resimap;
    old_resiptr = old_molptr->residues;
    old_count = old_molptr->res_count;
    while ( old_count-- )
    {
	if ( InResidueRange( old_resiptr ) )
	{
	    new_resiptr = new_molptr->residues;
	    new_count = new_molptr->res_count;
	    while ( new_count )
	    {
		if ( EQ( new_resiptr->res_id, old_resiptr->res_id ) )
		    break;
		new_resiptr++;
		new_count--;
	    }
	    if ( !new_count )
	    {
		ERROR( "Can't find residue id - Code error" );
		exit( 1 );
	    }
	    *mapptr++ = FindResidueIndexFromNum( new_resiptr->res_num_new, new_molptr );
	                                     /* new sequential residue index */
	}
	else
	    *mapptr++ = -1;
	old_resiptr++;
    }

/* print the data */

    if ( OptionTrue( opt_prtmap ) )
    {
	puts( "The atom map:" );
	mapptr = atommap;
	old_atomptr = old_molptr->atoms;
	old_count = old->atom_count;
	while ( old_count-- )
	{
	    WriteFAtom( stdout, old_atomptr++, use_both );
	    printf( ": %i\n", *mapptr++ );
	}
	puts( "The residue map:" );
	mapptr = resimap;
	old_resiptr = old_molptr->residues;
	old_count = old_molptr->res_count;
	while ( old_count-- )
	{
	    WriteFResidue( stdout, old_resiptr++ );
	    printf( ": %i\n", *mapptr++ );
	}
    }
    

/* set return data */

    *atom_map = atommap;
    *residue_map = resimap;
}
